library(dplyr) 
library(ggplot2) 
library(stargazer) 
library(tidyverse) 
library(gridExtra)

setwd("your_path_here")

ej <- read.csv("england_japan.csv", header=TRUE)

# Figure 1: 
# Percentage of respondents agreeing with 
# the three key statements across England and Japan.

jp1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$japan == 1]) / length(ej$manure1[!is.na(ej$manure1) & ej$japan == 1]) * 100, 2)
jp3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$japan == 1]) / length(ej$manure3[!is.na(ej$manure3) & ej$japan == 1]) * 100, 2)
jp5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$japan == 1]) / length(ej$manure5[!is.na(ej$manure5) & ej$japan == 1]) * 100, 2)
ep1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$japan == 0]) / length(ej$manure1[!is.na(ej$manure1) & ej$japan == 0]) * 100, 2)
ep3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$japan == 0]) / length(ej$manure3[!is.na(ej$manure3) & ej$japan == 0]) * 100, 2)
ep5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$japan == 0]) / length(ej$manure5[!is.na(ej$manure5) & ej$japan == 0]) * 100, 2)

q1 <- "Willing to eat food from HEBF..."
q3 <- "I have health concerns..."
q5 <- "Willing for public plants..."

Japan <- c(jp1, jp3, jp5) 
England <- c(ep1, ep3, ep5) 
Question <- c(q1, q3, q5) 

percentages <- data.frame(Question, Japan, England) 

# display those percentages in simple text format

percentages

# let's visualise that 

# convert the data from wide to long format
percentages_long <- tidyr::gather(percentages, Country, Percentage, -Question)

# convert 'Question' to a factor with the original order
percentages_long$Question <- factor(percentages_long$Question, levels = unique(percentages$Question))

# create the bar chart

ggplot(data = percentages_long, aes(x = Question, y = Percentage, fill = Country)) +
  geom_bar(stat = "identity", position = "dodge") +
  geom_text(aes(label = Percentage), vjust = -0.5, position = position_dodge(width = 0.9)) +
  facet_grid(~Country) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        plot.margin = unit(c(0.5, 0.5, 0.5, 3), "lines")) +
  labs(title = "Survey Responses by Country",
       x = "Question",
       y = "Percentage (%)")


# Table 1: Descriptive statistics

categories <- c(

	"EN Respondents", 
	"JP Respondents", 
	"EN Female (%)", 
	"EN Male (%)", 
	"JP Female (%)", 
	"JP Male (%)", 
	"EN Mean age", 
	"EN Max age", 
	"EN Min age", 
	"EN SD age", 
	"JP Mean age", 
	"JP Max age", 
	"JP Min age", 
	"JP SD age", 
	"EN Trust government (%)", 
	"JP Trust government (%)", 
	"EN Trust other people (%)", 
	"JP Trust other people (%)", 
	"EN Good health (%)",
	"JP Good health (%)",
	"EN University degree (%)",
	"JP University degree (%)"
) 

values <- c( 

	# EN respondents: 
	length(ej$id[ej$england==1]), 
	
	# JP respondents: 
	length(ej$id[ej$japan==1]), 
	
	# EN Female (%): 
	round(length(ej$id[ej$gender==2 & ej$england==1]) / length(ej$id[ej$england==1])* 100 , 2), 
	
	# "EN Male (%)",
	round(length(ej$id[ej$gender==1 & ej$england==1]) / length(ej$id[ej$england==1])* 100 , 2), 
	
	# "JP Female (%)", 
	round(length(ej$id[ej$gender==2 & ej$japan==1]) / length(ej$id[ej$japan==1])* 100 , 2), 
	
	# "JP Male (%)", 
	round(length(ej$id[ej$gender==1 & ej$japan==1]) / length(ej$id[ej$japan==1])* 100 , 2),
	
	# "EN Mean age", 
	round(mean(ej$age[ej$england==1]), 2), 
	
	# "EN Max age", 
	max(ej$age[ej$england==1]),
	
	# "EN Min age", 
	min(ej$age[ej$england==1]),
	
	# "EN SD age", 
	round(sd(ej$age[ej$england==1]), 2),
	
	# "JP Mean age", 
	round(mean(ej$age[ej$japan==1]), 2), 
	
	# "JP Max age", 
	max(ej$age[ej$japan==1]),
	
	# "JP Min age", 
	min(ej$age[ej$japan==1]),
	
	# "JP SD age", 
	round(sd(ej$age[ej$japan==1]), 2),
	
	# "EN Trust government (%)", 
	round(length(ej$trustGov[ej$trustGov > 4 & ej$england==1 & !is.na(ej$trustGov)]) / length(ej$trustGov[ej$england==1 & !is.na(ej$trustGov)])*100, 2), 
	
	# "JP Trust government (%)", 
	round(length(ej$trustGov[ej$trustGov > 4 & ej$japan==1 & !is.na(ej$trustGov)]) / length(ej$trustGov[ej$japan==1 & !is.na(ej$trustGov)])*100, 2),
	
	# "EN Trust other people (%)", 
	round(length(ej$trustOthers[ej$trustOthers > 4 & ej$england==1 & !is.na(ej$trustOthers)]) / length(ej$trustOthers[ej$england==1 & !is.na(ej$trustOthers)])*100, 2), 
	
	# "JP Trust other people (%)", 
	round(length(ej$trustOthers[ej$trustOthers > 4 & ej$japan==1 & !is.na(ej$trustOthers)]) / length(ej$trustOthers[ej$japan==1 & !is.na(ej$trustOthers)])*100, 2),
	
	# "EN Good health (%)",
	round(length(ej$ownHealth[ej$ownHealth > 4 & ej$england==1 & !is.na(ej$ownHealth)]) / length(ej$ownHealth[ej$england==1 & !is.na(ej$ownHealth)])*100, 2), 
	
	# "JP Good health (%)",
	round(length(ej$ownHealth[ej$ownHealth > 4 & ej$japan==1 & !is.na(ej$ownHealth)]) / length(ej$ownHealth[ej$japan==1 & !is.na(ej$ownHealth)])*100, 2),
	
	# "EN University degree (%)",
	round(length(ej$university[ej$university==1 & ej$england==1]) / length(ej$university[ej$england==1])*100, 2),
	
	# "JP University degree (%)"
	round(length(ej$university[ej$university==1 & ej$japan==1]) / length(ej$university[ej$japan==1])*100, 2)
	
) 

descriptives <- data.frame(categories, values) 

descriptives

# Figure 2: 
# Percentage of respondents agreeing 
# with the three key statements across 
# England and Japan, broken down by gender


jm1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$gender==1 & ej$japan==1]) / length(ej$manure1[!is.na(ej$manure1) & ej$gender==1 & ej$japan==1]) * 100, 2)
jm3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$gender==1 & ej$japan==1]) / length(ej$manure3[!is.na(ej$manure3) & ej$gender==1 & ej$japan==1]) * 100, 2)
jm5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$gender==1 & ej$japan==1]) / length(ej$manure5[!is.na(ej$manure5) & ej$gender==1 & ej$japan==1]) * 100, 2)
um1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$gender==1 & ej$england==1]) / length(ej$manure1[!is.na(ej$manure1) & ej$gender==1 & ej$england==1]) * 100, 2)
um3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$gender==1 & ej$england==1]) / length(ej$manure3[!is.na(ej$manure3) & ej$gender==1 & ej$england==1]) * 100, 2)
um5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$gender==1 & ej$england==1]) / length(ej$manure5[!is.na(ej$manure5) & ej$gender==1 & ej$england==1]) * 100, 2)

jf1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$gender==2 & ej$japan==1]) / length(ej$manure1[!is.na(ej$manure1) & ej$gender==2 & ej$japan==1]) * 100, 2)
jf3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$gender==2 & ej$japan==1]) / length(ej$manure3[!is.na(ej$manure3) & ej$gender==2 & ej$japan==1]) * 100, 2)
jf5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$gender==2 & ej$japan==1]) / length(ej$manure5[!is.na(ej$manure5) & ej$gender==2 & ej$japan==1]) * 100, 2)
uf1 <- round(length(ej$manure1[ej$manure1 > 4 & !is.na(ej$manure1) & ej$gender==2 & ej$england==1]) / length(ej$manure1[!is.na(ej$manure1) & ej$gender==2 & ej$england==1]) * 100, 2)
uf3 <- round(length(ej$manure3[ej$manure3 > 4 & !is.na(ej$manure3) & ej$gender==2 & ej$england==1]) / length(ej$manure3[!is.na(ej$manure3) & ej$gender==2 & ej$england==1]) * 100, 2)
uf5 <- round(length(ej$manure5[ej$manure5 > 4 & !is.na(ej$manure5) & ej$gender==2 & ej$england==1]) / length(ej$manure5[!is.na(ej$manure5) & ej$gender==2 & ej$england==1]) * 100, 2)

JPMen <- c(jm1, jm3, jm5) 
JPWomen <- c(jf1, jf3, jf5) 

EnMen <- c(um1, um3, um5) 
EnWomen <- c(uf1, uf3, uf5) 

pc2 <- data.frame(Question, JPMen, JPWomen, EnMen, EnWomen)


# Convert the data from wide to long format, separating by gender and country
pc2_long <- tidyr::pivot_longer(pc2, 
								cols = -Question, 
								names_to = c("Country", "Gender"),
								names_pattern = "(JP|En)(Men|Women)")


# Convert 'Question' to a factor with the original order
pc2_long$Question <- factor(pc2_long$Question, levels = unique(pc2$Question))

# Adjust Country names
pc2_long$Country <- recode(pc2_long$Country, JP = "Japan", En = "England")
	   
# Create bar chart with numbers 

ggplot(data = pc2_long, aes(x = Question, y = value, fill = Gender)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_text(aes(label = round(value, 1)), vjust = -0.5, position = position_dodge(width = 0.9)) +
  facet_grid(~Country) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
		plot.margin = unit(c(0.5, 0.5, 0.5, 3), "lines")) +
  labs(title = "Survey Responses by Gender and Country",
	   x = "Question",
	   y = "Percentage (%)")

# Table 2: 
# Support for HEBF across three different scenarios

# Flip manure3
# to turn it into 
# I do not have health concerns

ej$manure3 <- 8-ej$manure3


m1 <- lm(manure1 ~ women + japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m2 <- lm(manure1 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m3 <- lm(manure3 ~ women + japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m4 <- lm(manure3 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m5 <- lm(manure5 ~ women + japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m6 <- lm(manure5 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)


# Display the output of the models
stargazer(m1, m2, m3, m4, m5, m6, type="text")

# Supplemental 
# Not reported in the paper. 
# Because there is a gender bias in UK data, 
# let's run some UK models with and without weighting

# create england only data frame 

en <- ej[ej$england == 1, ]

m1 <- lm(manure1 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en)
m3 <- lm(manure3 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en)
m5 <- lm(manure5 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en)

m2 <- lm(manure1 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en, weights = W8)
m4 <- lm(manure3 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en, weights = W8)
m6 <- lm(manure5 ~ women + age + trustGov + trustOthers + ownHealth + university, data = en, weights = W8)

stargazer(m1, m2, m3, m4, m5, m6, type="text")

# Figure 3: 
# Interactions

# Rebuild models on Japan and England data 

ej$women <- factor(ej$women) 
ej$japan <- factor(ej$japan) 

m2 <- lm(manure1 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m4 <- lm(manure3 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)
m6 <- lm(manure5 ~ women * japan + age + trustGov + trustOthers + ownHealth + university, data = ej)




predict_margins <- function(model, data) {
  # Create a new data frame for prediction that contains all combinations of 'women' and 'japan'
  new_data <- expand.grid(women = levels(data$women), japan = levels(data$japan))
  
  # Add the mean values for other variables
  new_data$age <- mean(data$age, na.rm = TRUE)
  new_data$trustGov <- mean(data$trustGov, na.rm = TRUE)
  new_data$trustOthers <- mean(data$trustOthers, na.rm = TRUE)
  new_data$ownHealth <- mean(data$ownHealth, na.rm = TRUE)
  new_data$university <- mean(data$university, na.rm = TRUE)
  
  # Calculate predicted values
  predictions <- predict(model, newdata = new_data, type = "response", se.fit = TRUE)
  new_data$predicted <- predictions$fit
  new_data$se <- predictions$se.fit
  return(new_data)
}


# Generate the plot for given predicted margins
plot_margins <- function(predicted_margins, title) {
  ggplot(predicted_margins, aes(x = women, y = predicted, color = japan, group = japan)) +
    geom_line(position = position_dodge(width = 0.1)) +
    geom_point(position = position_dodge(width = 0.1)) +
    geom_errorbar(aes(ymin = predicted - 1.96 * se, ymax = predicted + 1.96 * se), 
                  width = 0.2, position = position_dodge(width = 0.1)) +
    scale_color_manual(values = c("0" = "blue", "1" = "red")) +  # Assign colors
    labs(title = title, x = "Women", y = "Linear Prediction") +
    theme_minimal()
}

margins_m2 <- predict_margins(m2, ej)
margins_m4 <- predict_margins(m4, ej)
margins_m6 <- predict_margins(m6, ej)

# Create the plots
p2 <- plot_margins(margins_m2, "Eat food")
p4 <- plot_margins(margins_m4, "No health concerns")
p6 <- plot_margins(margins_m6, "Use in public parks")

combined_plot <- grid.arrange(p2, p4, p6, nrow = 1)



